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AIRBORNE RADAR TECHNOLOGY FOR UINDSHEAR DETECTION 


By 

Joseph L. Hibey* and Camille S. Khalaf** 

ABSTRACT 

The objectives and accomplishments of the two and a half year effort to 
describe how returns from an on board Doppler radar are to be used to detect 
the presence of a wind shear are reported. The problem is modeled as one of 
first passage in terms of state variables, the state estimates are generated 
by a bank of extended Kalman filters working in parallel, and the decision 
strategy involves the use of a voting algorithm for a series of likelihood 
ratio tests. The performance issue for filtering is addressed in terms of 
error-covariance reduction and filter divergence, and the performance issue 
for detection is addressed in terms of using a probability measure trans- 
formation to derive theoretical expressions for the error probabilities of a 
false alarm and a miss. 


♦Associate Professor, Department of Electrical and Computer Engineering, Old 
Dominion University, Norfolk, Virginia 23529. 

♦♦Graduate Research Assistant, Department of Electrical and Computer Engi- 
neering, Old Dominion University, Norfolk, Virginia 23529. 
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1. INTRODUCTION 


We shall begin by briefly summarizing the goals and accomplishments of 
the two and a half year duration of this grant. Following this overview, we 
shall present a more detailed exposition and then conclude with some ideas 
for future work. 

The primary objective involved detecting the presence of wind shear and 
alerting the pilot in time to avoid catastrophe. To accomplish this, we 
have assumed throughout that a Doppler radar would reside on board the air- 
craft and provide us with measurements, or returns, that we could then pro- 
cess for purposes of detection and ultimate avoidance. The feasibility of 
such an approach was indicated in a 1983 report [1] of the National Research 
Council on wind shear. 

With this in mind, work proceeded by first attempting to obtain a 
mathematical description of the problem. This was achieved by representing 
the onset of wind shear as a first-passage problem. Also called the 
classical disruption problem or exit problem, it essentially involves 
monitoring the evolution of a stochastic process and attempting to determine 
the first time it exceeds some threshold. The actual models were described 
in terms of state variables that satisfied a system of stochastic 
differential equations. The radar returns were represented as nonlinear 
signals in additive white Gaussian noise. Nonlinear filtering techniques of 
the extended Kalman type were then applied to estimate both the time of 
occurrence of wind shear and the magnitude of the so-called microburst. 

This theoretical formulation and ideas for its simulation and eventual im- 
plementation formed the core of the first year and half effort, and is fully 
described in Ref. [2]; it also served as a progress report for that period 
and is included in Appendix A to this document. 



Encouraged by the preliminary findings of the initial research, we 
continued in the last year by refining our models to more appropriately 
reflect reality and by addressing the issue of performance. Thus, we first 
explicitly took into account phase information and modified the filtering 
algorithms accordingly. Many computer simulations were run to determine 
viable approximations leading to acceptable results. We then shifted our 
attention to the system's overall performance in terms of evaluating the 
false alarm and miss error probabilities. Although the findings are of a 
theoretical nature, we shall include some ideas regarding their integration 
into an actual algorithm. 

In the next section we shall formulate the problem in terms of a math- 
ematical model involving state space representations and use nonlinear fil- 
tering theory to derive algorithms for estimating the quantitites needed to 
detect wind shear. Computational issues related to simulation and implemen- 
tation procedures will then be discussed. Finally, this will lead us to a 
presentation of detection and performance issues to be followed by some 
concluding remarks. 

Before proceeding, we mention that in the following presentation we 
shall often refer to paper [2] in Appendix A for notation and some of the 
basic techniques we have adopted. As will be seen, in some instances the 
equations used there will not be modified, while in other cases they are. 
Therefore, for the sake of completeness and ease in reading, if an equation 
has not changed, it will merely be repeated in this report with the same 
numbering system and with little or no comment. On the other hand, variants 
of the original equations will appear in this report with a "prime" synbol 
attached to the equation number. Finally, starting with Eq. (12) in this 
report, there are no such corresponding equations in the earlier report. 
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2. MATHEMATICAL REPRESENTATIONS 


Modeling 

In [2], we concluded that the radar return can be modeled as 


r t = h (s t> + v t 

where h(sj = A sin {a> + 4ir s / xl t, 

t 1 c t 

s t = b t z t + x t (1 ‘ z t } 


( 1 ) 


( 2 ) 


and 


z t 



if t < T 
if t > T 


(3) 


This is the classic "signal h (• ) in additive white Gaussian noise v" prob- 
lem. It is somewhat unconventional, however, in that it involves z^, a 
zero-one process that distinguishes between the absence of wind shear (z = 0 
for t < T, implying a relatively low wind velocity s = x) and the presence 
of wind shear (z = 1 for t > T, implying a relatively high wind velocity 
s = b). The idea of detection, then, is to derive an estimate T of T, the 
first time of occurrence of a microburst. As such an estimate will require 
an estimate z of z, one must first find state-space models for all the 
processes involved and then derive an appropriate filter. 

The earlier model describing the states has now been modified to 
include phase information <|». As usual, we shall assume that i|> is a 
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uniform random variable distributed over (-*, ir). The resulting model then 
becomes 
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observation: dy t = A sin {oj c t + 4 * ^ [b t z t + x t ( 1 -z t dt + dV t (5') 

v .. «, 


where U. is Brownian motion (mean zero and variance U t) that is mutually 
t c 

independent of the above disturbances and initial conditions; with ^ chosen 
as uniform, H and U will typically be selected very small. (Note: in [2], 
the wavelength v in (5‘) was ambiguously written as a). 

Signal Statistics 

This new measurement equation involving the process ^ can be given the 
following interpretation. In actuality, at a given range gate there exists 
a large number of scatterers each with a different velocity. If the sam- 
pling rate is sufficiently high, then each of these velocities is approxi- 
mately constant. Therefore, one can model this situation by imagining one 
scatterer at each range gate having a velocity that varies with sample time. 
What is important is that the model possess the same statistics as the real 
system, and this is accomplished by requiring a uniform distribution for the 
random phase process Given this model, the filters used to derive tne 
appropriate state estimates will now include additional equations to esti- 
mate as well . 

The signal h(») can be viewed as an example of exponential modulation 
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(see, for example. Gray and Davisson [3, pp. 258-260] or Weinstein and 
Zubakov [4, Chap. 8)), that is, 


4irt 


h(t) = Im q(t) = Im (A exp [j(a> c t + [b t z t + x t (l-z t ) + i|>)]} 


where Im {•} denotes the imaginary part of {•}. From this, the statistics 
of h can be derived from those of q. Thus, since the amplitude A has a 
Rayleigh distribution, the velocity x (or b) has a Gaussian distribution, 
and the phase i|> has a uniform distribution, it follows that, given z (0 or 
1), the sinusoid in h is distributed as 1 /tt ^1 - y^ for | y | <1 (see Papoulis 
[5, p.100]) and h is Gaussian (see Doviak and Zrnic [6, p. 50]). The 
results of Doviak and Zrnic [6, p. 440] and Papoulis [5, pp. 268, 269] then 
lead to an autocorrelation R of q of the form 

q 


V* > ■ 


c e 


■dr 


J V 


where c and d are constants, and this implies its power spectrum is 
Gaussian (also see Doviak and Zrnic [6, p. 43]). 

The previous result is valid conditionally, that is, with constants c, 
d corresponding to x^ when z t = 0 (t < T) and there is no wind shear, and 
with constants c' , d' corresponding to when z^ = 1 (t>T) and there is 
wind shear. To derive an expression valid for all t, one finds through 
reconditioning and the fact that T is exponentially distributed (see [2] in 
Appendix A) that R^ is of the form 
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V T > - 


-At 


E e 


ja(x t - 


X ) 
s' 


fl-e“ Xt ) 


E e 


ja' (b t -b $ ) 


( 6 ) 


with a and a' constants. But since the increments x. - x and b, - b are 

t s t s 

themselves Gaussian, we find as above that 


R q( T ) = 


-At 


-dx 


c e 


JV 


(l-e' U ) c 1 



J'V 

e 


(7) 


which shows that q, and therefore h, are nonstationary processes. This 
result suggests the use of Kalman-type filtering which, unlike Wiener 
filtering, makes no assumptions regarding the necessity of stationarity. 
Filtering 

Using this model, we next generate MMSE estimates (denoted with a ~ 
symbol) of the state variables and use them to detect the onset of wind 
shear. Thus, as in the earlier report [2], we find that the conditional 
mean z t given by 

z t = Prft > T | y $ , 0 < s < t] (8) 


represents the probability that wind shear has occurred. But again, 
recognizing the need for a suboptimal estimate z^ because of the nonlinear 
nature of the problem, a suboptimal estimate T of the microburst occurrence 
time T is given by 


T = min { t j z t > k} 


(9) 
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where k is some threshold in the unit interval. This characterizes our 
problem as one of first-passage. 

Our next step is to specify the filtering algorithms required to 
generate the estimates used in detection. These will be based on the 
extended Kalman filter and, for the system in (4') and (5‘) above, are given 
as follows: 



where IN = dy - (h + — . SOT) dt. Innovations 

2 

A 

SOT = trace [PH], Second-Order-Term, and 
A 2 - 2 

H = 3 h / 3 (T) , Hessian matrix. 
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Until now, all we have done is modify the appropriate equations in [2] to 
account for phase information. Looking ahead, however, to our discussion 
about performance, it will be helpful to explicitly represent the gain 

_ _ A 
P 9 h / a (• ) in (10' ) and h above, with <f> = 4tt / v , as 


9h/9x G 

A x A 

P a h/3 z = G z 

ah/ab G, 

b 

a h/a^T g^ _ 


Pll 

p 12 

P 13 

p 14 

P 21 

P 22 

P 23 

p 24 

p 31 

P 32 

P 33 

P 34 

P 41 

p 42 

p 43 

p 44 


b - z <j> t A cose 
z 

( K )' 1 


where 


A 

h = A sin {u> c t + <j>t [b t z t + x^ (1-z^.)] + ^ } = A sin e . (13) 

Rewriting (10'), we therefore obtain the filter 

dx = -F 7 dt + G x V c _1 - IN 

dz = X (1-7) dt + G z V' 1 . IN 

db = -G b dt + G b V" 1 - IN 

dJT = -H f dt + fy V" 1 . IN 

3. COMPUTATIONAL ISSUES 

The above formulation has been given in germs of continuous time. A 
discrete-time formulation was derived from this using standard Euler 
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approximations (see [2]) for purposes of simulation on a digital computer. 

As was expected, because of the nonlinear nature of the measurement equation 
given in (5 1 ), the estimates produced by the filters in (10 1 ) sometimes 
diverged for certain choices of parameter such as initial conditions 
(x , b , P , etc.), noise intensities (W , B , V , etc.), and signal power 

U U U v v v 

levels ( related to A in (5')). More precisely, the error-covariance matrix 
P in (ll 1 ) did not remain nonnegative definite as it should. The reasons 
why this occurs include the approximation techniques utilized in determining 
an extended Kalman filter and the finite word length or limited precision of 
a digital computer; in the latter case, truncation errors tend to accumulate 
over time (see Jazwinski [7, chap. 9]). 

In order to rectify this situation, the following approaches were 
tried. The diagonal elements of the P matrix in (IT) represent error- 
variances of the appropriate state estimates. Being variances, these should 
never become negative. Therefore, as a first attempt, these variables were 
monitored and reset to zero each time they became negative. 

With only marginal improvement from the above, we considered the off- 
diagonal elements of P which represent cross error-covariances. For 
example, the typical (i,j) element for i * j is of the form 






where E denotes expectation given the measurements, 


inequality, an upper bound on p . results in 

J 


Using the Schwarz 


ij 


[E(u n 


/N \2 

■ u i) 


E <yV 2 ] 


1/2 


= fPil 


P jjJ 


1/2 


(15) 
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Thus, by monitoring the p..'s, we were able to reset values to insure the 

1 J 

validity of (15) for all time. 

Although this approach showed further improvement, it was still not 
satisfactory. At best, it is an ad hoc technique that works in only limited 
cases. We therefore turned to square-root filtering algorithms (see, for 
example, Anderson and Moore [8, chap. 5]). These algorithms are based on 
so-called singul ar-value-decompositions and are specifically designed to 
maintain the nonnegative definite property of P and thereby eliminate filter 
divergence. Although their use implies a more complex filter, the improve- 
ment in the quality of the estimates justifies their adoption. A copy of a 
computer program that generates state estimates from such a square-root 
filter of the Kalman type appears in Appendix B with some examples of its 
output . 

Turning now to the implementation aspect of the problem, we refer once 
again to [2] for the basic methodology. There we see the need for a bank of 
N of the above type filters each receiving its sample measurements from a 
particular range location, or range gate, within the region affected by the 
microburst; each filter will process a batch of M measurements in a serial 
fashion (the results in Appendix B use M=32), and all filters will work in 
parallel. The actual decision as to whether or not wind shear has occurred 
will be based on the outputs z of these N filters, and our approach will 
employ a voting algorithm to be explained in the next section on detection 
and performance. 

4. DETECTION AND PERFORMANCE 

The performance issue involves determining how good the filtering and 
detection algorithms work. In the case of filtering, this is measured by 
examining the errors in estimation and is directly reflected in the elements 
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of the error-covariance matrix P in (IT). Whereas this was the subject of 
the previous section, we now address performance in the context of 
detection. 

A measure of how good a detection scheme works involves an assessment 
of the errors one commits while making decisions. Here we are confronted 
with two types of errors: declaring wind shear when it is not present, and 
declaring no wind shear when indeed it is present. These errors, called 
false alarm and miss, respectively, are often very difficult to compute. 

Our objective in this section, then, is to explain how one might evaluate 
them both theoretically and empirically. 

Likelihood-Ratio Test 

Our strategy for deciding about the presence or absence of wind shear 
is based on a likelihood ratio test such as described by Van Trees [9, chap. 
2]. As indicated above, we are interested in determining the time at which 
the process z t in (3) changes abruptly from zero to one; the detection of 
such abrupt changes is discussed in the text edited by Basseville and 
Beneveniste [10]. As it turns out, the likelihood ratio involves the 
estimate z^ of z^ given in (8), or, more precisely, its suboptimal estimate 
z t given in (14). We therefore define the likelihood ratio L as 

A Pr [wind shear | observations of y from 0 to t] 

L = 

Pr [no wind shear | observations of y from 0 to t] 

and use (8) with z replaced by *z” to obtain 
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z 


(16) 


A Pr [T < t j y , 0 < s < t] 

L = ! 

Pr [T > t | y s> 0 < s < t] 1-7 


The likelihood-ratio-test (LRT) involves testing the magnitude of L in 
relation to some threshold y' defined as 


y = 


A p n ( C 10 _( W 


P 1 < c or c n> 


(17) 


where Pq and P^ are the a priori probabilities of the absence and presence 
of wind shear, respectively, and C.., i ,j = 0,1 are the relative costs of 

’ J 

making correct (i=j) and incorrect (i*j) decisions. For example, as is 
often the case, if no cost is attached to making correct decisions, C Q0 = 

= 0; a choice of C 1Q = 3, say, implies a cost of deciding in favor of wind 
shear (i=l) when indeed there is none (j=0), while a cost of C Q1 = 19, say, 
is incurred in deciding against wind shear (i=0) when in fact it is present 
(j=l). Also, by way of example, Pq = .95 indicates a strong belief that no 
wind shear is present and therefore implies P^ = .05. The actual values 
used in any test will be chosen by the pilot based on airline safety 
procedures and standards, and current weather reports. From the choices 
selected above, the interpretation is that one is less tolerant of missing a 
wind shear (Cqj> C^q) in a situation where wind shear appears to be unlikely 
(Pq > > P^). We therefore arrive at the following LRT: 
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LRT : 


> 

< 

H, 


(18) 


Here, H Q denotes the hypothesis there is no wind shear and denotes the 
hypothesis there is wind shear. Equivalently, assuming y ' > 0, we have 

H i 

LRT: z < y ( 19) 

H 0 

A - 

where y = y' / (1+y 1 ). Thus, we compute z and declare the presence of wind 
shear if z > y, and declare its absence otherwise. 

Voting Algorithm 

We next recall that in reality we have a bank of N filters each using 
measurements taken at a specific range gate. Therefore, we will need N 
LRT's and some type of voting algorithm on which to base our decision. We 
proceed as follows. 

Let us assune a synmetric microburst pattern and, for convenience, N + 
1 filters acting over a range that encloses the resolution volume of the 
weather target, as shown below: 



A 

Here, range time t s = n At, n = - N/2, ..., -1, 0, 1, .... N/2; also see 
Fig. 3 in [2]. We further assume that if, say, 4/5 of the tests indicate 
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wind shear (H^), we declare wind shear; if 4/5 of the tests indicate no wind 

shear (Hq), we declare no wind shear; otherwise, we make no decision at this 

time and wait for more data. To express this, we first denote the indicator 

function If.) of an event A by 
( 1 if A occurs 


I I A] = 


J 0 otherwise 

We next define the ninber J of "triggered gates" (those indicating wind 

shear) by 

A N/2 

J = l I[z t (iAx ) > y] (20) 

1 =-N/2 X 


_ A 

and its complement by J = (N+l-J) (the number of gates indicating no wind 
shear). Then, setting the ratio r = 4/5, our voting algorithm becomes 


if J > r (N+l), then decide H^; 
if "J > r (N+l), then decide Hq; 
otherwise, make no decision. 


( 21 ) 


The ratio r used in (21) is a design parameter set by the pilot in 
accordance with weather reports and existing standards and guidelines. For 
example, with a weather report indicating virtually no winds, 4/5 might be 
viable because any microburst would probably be symmetric; on the other 
hand, a report of a 20 knot wind might suggest the use of ratio 2/3 because 
any microburst would probably be asynmetric. In general, as weather reports 
indicate higher wind velocities, a lower number of "gates" need to be 
"triggered" before declaring the presence of wind shear because its 
occurrence becomes more probable. Finally, we note that z is a process 
evolving with time t, so any meaningful test would actually require that 
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(21) give consistent results over some prespecified duration of time before 
the ultimate decision is made; examining over this duration would also 
provide a direct correlation with this decision. 

Error Analysis 

Referring to (19), we see that a so-called type 1 error or false alarm 

is committed if hypothesis Hq (no wind shear) is indeed true, yet we 

mistakenly conclude that hypothesis (wind shear) is true because z > y. 

Likewise a so-called type 2 error or miss is committed if is true, yet we 

decide in favor of Hq because z < y. Therefore, denoting the conditional 

densities of z by p(z|H. ), i=0, 1, we define the error probabilities of a 

false alarm P.. and miss P,. by 
r M 

A 1 

P F = Pr [z t > y I Hq] = / p(z t I Hq) d z t (22) 

Y 

A _ Y _ _ 

p M = p r [ z t < Y | Hj = / p (z t | H x ) d z t . (23) 

One approach to computing Pp and P^ is to use Monte Carlo techniques. 
For example, considering Pp (similar remarks apply to P M ), one would gener- 
ate observations from (5‘), but with z^ = 0 which corresponds to Hq being 
true as in (22). The observations would then be used in (10') and (IT) to 
evaluate z t , and finally Pp could be evaluated using relative frequency type 
calculations. 

On the other hand, (22) and (23) show the need to evaluate the condi- 
tional densities p (z j H^), i=0, 1 of the likelihood ratio z, but quite 
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often this is not possible. In such cases we might be satisfied with ob- 
taining, for example, upper bounds on these error probabilities as in [11]. 
Although we shall not follow that approach here, some of the ideas in [11] 
will be used to derive theoretical expressions for and P^. If they could 
then be evaluated, their values could be compared with the empirical values 
obtained from the Monte Carlo techniques mentioned above. We now proceed 

with our analysis of P.. defined in (22); similar comments apply to P.. in 

r M 

(23) and, for most part, will be omitted. 

Our approach can be summarized as follows. If of (10') were an 
optimal estimate of z^. rather than a suboptimal estimate, then the 
innovations process IN driving (10') would be a Brownian motion process and 
the density p (z) of z t would satisfy the Fokker-Pl anck (FP) equation 

— p(z)] + [ (G V' 1 ) 2 p(z)l ; ’ (24) 

— ? 1 L 

3 1 3 z 2 3z 

this is also called the Kolmogorov forward equation. However, this equation 
does not exist because z t is in fact a suboptimal estimate of z^, so IN is 
not Brownian motion. Furthermore, (22) requires the conditional density 
p(z | Hq), not the unconditional density satisfying (24) which has 
coefficients A(l-z) and G z V c * that are the same regardless of which 
hypothesis is true. For these reasons, therefore, we shall use a measure 
transformation (see [11] and [12] for a similar application) which will lead 
to a FP equation with coefficients that do indeed distinguish between the 
two hypotheses. The cost, however, will be an increase in the dimension of 
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the states whose density satisfies this equation, so the conditional density 
required in (22) will follow by integration. 


To begin the derivation, we define the martingale 
a t _ 

m. = / f(h) [dy - h dr] (25) 

u o t T 

where f(h) will be specified shortly. Then defining the exponential formula 
(see van Schuppen [13, pp. 32, 33] or Brtmaud [14, pp. 336-338]) 

A t „ 

e(m.) = exp [m t - / f (TT) V dr] 
c t Q c 

t _ t 9 

= exp {/ f(h) [dy - h d t] -/ r(F) V d x} (26) 

0 0 c 


it follows that if E e (m^) = 1, then e(m t ) can be used with the original 
probability measure Pq corresponding to Hq to define a new probability 
measure (see van Schuppen [13, p. 38] or Lipster and Shiryayev [15, pp. 
314, 315]) such that 

e(m t ) = dP 2 /dP Q . (27) 

The verification that E e(m T ) = 1 follows from van Schuppen [13, pp. 34-36] 

or Lipster and Shiryayev [15, p. 295] by invoking a sufficient condition 

2 -1 
requiring that f (TT) be bounded; thus choosing f(TT) = 7iV c satisfies the 

requirement because h is a sinusoid bounded by 1. Turning once again to van 
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Schuppen [13, pp. 40-42], we use the martingale translation theorem to 
obtain a measurement equation with respect to the new measure given by 

dy = h dt + h dt + dM (28) 

where the martingale H is a Brownian motion process; for a proof of this 
last assertion, see [16, p. 53] for a similar derivation. Finally, we use 
(28) in (14) to obtain the new filtering equations 

d7 = [ -F 7 + G x V^ 1 h] dt + G x V^ 1 dM 

dz = [A (1-7) + G z V' 1 h] dt + G z V^ 1 dM 

db = [ -G¥ + Gjj V" 1 h] dt + G fa V^ 1 dM 

dT = [ -H f + G^ V” 1 h] dt + G^ V" 1 dM 

As (29) indicates, the "drift" coefficients now explicitly contain the 
signal h of (13) and, since h depends explicitly on z, these coefficients 
can distinguish between the hypotheses and as z goes from 0 to 1, 
respectively. However, because our interest is currently focused on 
evaluating Pp, we set z = 0 and express h by 

A 

h = A sin (u> c t + x^ + = A sin e Q . (30) 

Therefore, substituting (30) into (29) indicates that the filter in (29) 
represents a coupling of not only the state estimates x, z, b, and but 
also of the original (unfiltered) states x of (4-a* ) and ip of (4-d 1 ). This 
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then is the anticipated increase in dimension alluded to above and leads to 
the following augmented system of equations: 


dx = [ -F7 + G x V " 1 h] dt + G x V ” 1 

dz = [A ( 1 - 7 ) + G z V ' 1 h] dt + G z V ^ 1 

db = [ -Gb + G b V ' 1 h] dt + G b V " 1 

cty = [-H f + G, V ' 1 h] dt + G, V " 1 

L ip c J ip c 

dx = -F x dt + d W 

dip = -H <j> dt + d U 


dM 

dM 

dM 

dM 


(31) 


This can also be expressed in vector form by defining the augmented state as 

A j 

x = [x z b i(i x i|»] and writing 


dx = A(x, t) dt + B (x, t) dW 


(32) 


where the column vectors A^t), B(x_,t) and W_ are defined accordingly. 

We now recall that our objective was to obtain an expression for 
p (z^ | Hq ) as required by (22). Therefore, using h of (30) in (32), the 
conditional joint density of _x given Hq, denoted by p(_x, t | Hq), satisfies 
the Fokker-Pl anck equation 


IP = _ 7 . (Ap) + i l 

at - 2 i,j=l 


3 2 ( V> 


3 x, 3 x. 

J 


(33) 


where v x denotes the gradient operator with respect to the vector x_ and 
^ T 

b = (B B ) ^. The fact that (33) is valid follows from our use of the 

J * J 
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measure transformation that resulted in the state equation (32) being driven 
by a Brownian motion process W_. The density p(z t | Hq) now follows from 

integrating p(_x, t | Hq), that is, 

00 00 IT __ 00 W 

p(z t I h q) = f dx f db f <u > f dx f d * p (x (t) I Hq) , (34) 

-® -» -ir -» -tr 

and the final expression for the probability of a false alarm defined in 
(22) becomes 


1 

P F = / d z f p (x, t | H q ) d(x, b, v, x, <(0 (35) 

y 

where the latter integration symbol denotes the multiple integration shown 
explicitly in (34). 

A similar analysis leads to an expression for the probability of a miss 
P M> To derive it, we set z = 1 in (13) and express h by 

& 

h = A sin { to t + 4> t b, + ip . > = A sin 0. . (36) 

' c t t 1 

Its use in (29) then represents a filter with coupling among the states x, 
z, b, and ip , and also the states b of (4-c* ) and i|> of (4-d 1 ). Therefore, 
the equation for x in (31) is replaced by the equation for b and the 
augmented state vector becomes jx = [x z b b The vectors A and B of 

(32) change accordingly, and (34) is replaced by 
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00 


00 


p(z t I H 1 ) = / d X f db r d ip f db f ckp p(_x, t | Hj). (37) 

-00 -« -IT -00 “IT A 


Then the final expression for the probability of a miss defined in (23) 
becomes 


Y _ _ _ _ 

P M = / d z / p(_x, t | H.) d(x, b, ip , b, ij>) (38) 

m 0 1 

To surrmarize what we have done, the false alarm error probability P 
and miss error probability P u each require the solution of the Fokker-Planck 
equation of (33) followed by the integrations in (35) and (38), 
respectively. Since an exact solution is virtually impossible, some 
nunerical procedure must be used to obtain suitable approximations. 
Considering a worst case scenario, however, by assuming large measurement 
noise such that ■+ 0), causes many terms in A and B to vanish and thereby 
makes (33) less formidable. 


5. CONCLUSIONS 

From the outset, we proposed to describe how one might process the 
returns from an on board Doppler radar system so as to detect the presence 
of a wind shear. Questions relating to specific hardware such as antenna 
design and so on were not the major concern. In contrast, the scope of tnis 
effort was to address software issues in the context of deriving filtering 
and detection algorithms to achieve the above objective. 

That this is a final report should in no way suggest the completion of 
this phase of the overall effort. A good part of the preceding formulation 
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has resulted in theoretical type expressions which ideally could be used as 
a basis for practical implementation, and we have endeavored to indicate 
this in the appropriate places of this report. It is our sincerest hope 
that the adoption of some of these ideas will be useful in the attempt to 
further understand and ultimately solve the wind shear problem. 

ACKNOWLEDGEMENT 

This is a final report on the research project "Airborne Radar 
Technology for Windshear Detection." Specific efforts were directed in the 
area of "Wind Shear Modeling and Detection: An Analysis in Terms of a 
First-Passage Problem." The performance period for this research ended 
August 31, 1988. This work was supported by the NASA Langley Research 
Center through NASA Research Grant NAG-1-626 and was monitored by Dr. Leo D. 
Staton, Technical Monitor of the GCD-Antenna & Microwave Research Branch, 
NASA Langley Research Center, MS/490. 


22 



REFERENCES 


1. National Research Council, Report of the Committee on Low-Attitude 
Wind Shear and its Hazard to Aviation, National Academy Press, 
Washington, D.C., 1983. 

2. C. S. Khalaf, J. L. Hibey, and L. D. Staton, "Wind Shear Modeling and 

Detection: An Analysis in Terms of a First-Passage Problem, presented 

at the Twentieth Asilomar Conference on Signals, Systems, and 
Computers, Pacific Grove, CA., Nov. 1986. 

3. R. M. Gray and L. D. Davisson, Random Processes : A Mathematical 
Approach for Engineers , Prentice-Hall, 1986. 

4. L. A. Wainstein and V. D. Zubakov, Extraction of Signals from Noise , 
Prentice-Hall, 1962. 

5. A. Papoulis, Probability, Random Variables and Stochastic Processes , 
McGraw-Hill, 1984. 

6. R. J. Doviak and D. S. Zrnic, Doppler Radar and Weather Observations , 
Academic Press, 1984. 

7. A. H. Jazwinski, Stochastic Processes and Filtering Theory , Academic 
Press, 1970. 

8. B. D. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, 

1979. 

9. H. L. Van Trees, Detection, Estimation, and Modulation Theory, Part I, 
Wiley, 1968. 

10. M. Basseville and A. Benvemste (ed.). Detection of Abrupt Changes in 
Signals and Dynamical Systems , Springer-Verlag, 1980. 

11. J. L. Hibey, D. L. Snyder, and J. H. van Schuppen, "Error Probability 
Bounds for Continuous-Time Decision Problems," IEEE Trans, on 
Information Theory, Vol. IT-24, No. 5, Sept. 1978, pp. 608-622. 

12. J. L. Hibey, "Cycle Slipping in an Optical Communication System 
Employing Subcarrier Angle Modulation," IEEE Trans, on Information 
Theory, Vol. IT-33, No. 2, Mar. 1987, pp. 203-209. 

13. J. H. van Schuppen, "Estimation Theory for Continuous Time Processes, 
a Martingale Approach," Ph.D. Dissertation, Umv. of Calif., Berkeley, 
1973. 

14. P. BrSmaud, Point Processes and Queues: Martingale Dynamics, 

Springer-Verlag, "1981. 

15. R. S. Lipster and A. N. Shiryayev, Statistics of Random Processes II , 
Applications , Springer-Verlag, 1978. 


23 



16. J. L. Hi bey, "Error-Probability Bounds for Continuous-Time Decision 
Processes," D.Sc. Dissertation, Washington University, St. Louis, 
1976. 


24 



APPENDIX A 


DEPARTMENT OF ELECTRICAL & COMPUTER ENGINEERING 
COLLEGE OF ENGINEERING & TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23508 


WIND SHEAR MODELING AND DETECTION: AN 

ANALYSIS IN TERMS OF A FIRST-PASSAGE PROBLEM 


By 

C. C. Khalaf 
and 

J. L. Hibey, Principal Investigator 


Progress Report 

For the period ending November 30, 1986 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

NASA Grant NAG-626 
Dr. Leo Staton 

FED-Antenna & Microwave Research Branch 


Submitted by the 

Old Dominion University Research Foundation 
P. 0. Box 6369 
Norfolk, VA 23508 


January 1987 



WIND SHEAR MODELING AND DETECTION: 

AN ANALYSIS IN TERMS OF A FIRST-PASSAGE PROBLEM 

By 

C. S. Khalaf 1 and J. L. Hibey 2 
SUMMARY 

The enclosed manuscript was presented at the Proceedings of the 1986 
Asilomar Conference on Signals, Systems, and Computers, which was held in 
Pacific Grove, California, in November 1986. 

The project work is supported by NASA Grant NAG-1-626. 


graduate Research Assistant, Department of Electrical & Computer 
Engineering, Old Dominion University Research Foundation, Norfolk, VA 
23508. 

2 Associate Professor, Department of Electrical & Computer Engineering, Old 
Dominion University Research Foundation, Norfolk, VA 23508. 


26 



1 


WIND SHEAR MODELING AND DETECTION: AN ANALYSIS IN TERMS OF A FIRST-PASSAGE PROBLEM 


* ★ ★★ 

C. S. Khalaf , J. L. Hibey and L. D. Staton 

* 

Old Dominion University, Norfolk, Virginia 
MS 490 NASA Langley Research Center, Hampton, Virginia 


«dLJ?A_CJ 

An aporcacn is proposed for extracting Infor- 
mation from airborne Ouppler radar returns that 
:an be used to determine the presence of wind 
shear in commercial aviation. Issues that are 
discussed nclude modeling in terms of stochas- 
tic di ffei-ential equations, estimation In terms 
}r minimummean-square-error filters, and detec- 
tion in terms of a first passage problem. Some 
discussion of the hardware architecture needed 
i process the data in a timely and efficient 
tanner 15 also included 


1 1 NTRQQUCT ION 

The problem of wind shear in commercial avia- 
tion has always existed, but on I v recently has it 
been given the attention it deserves In the last 
f’l-teen vears. it has been listed as the cause of 
sc eial accidents involving large carriers, with 
*-ne most recent being the crash of Delta Airlines 
c . 1 ght 191 in Dallas on August 2, 1985 that 
claimed the lives of 137 people 

The phenomenon can be loosely described as a 
jU'id-’n downburst of wind at a velocity of aproxi- 
matel 40 meters/second which, upon impact with 
the urcund. results in a hor'zontal component of 
wirq velocity with high magnitude, and with 
diiecrion that aorupt 1 y changes as one moves 
through it at an altitude close to the ground, 
thus it poses a severe hazard to aircraft who 
encounter this so-called microburst (see fujita 
(II) in either the take-off or landing phase of 
flight. In that the aircraft first experiences a 
high head wind causing lift followed suddenly by a 
high tail wind causing descent. Therefore, with 
little altitude in which to maneuver, pilots are 
unable to recover control and disaster results 

The manv causes of wind shear are discussed in 
report [2], also see Doviak and Zrnlc (3, chap 3] 
These range from wet microbursts that occur in the 
wake of thunderstorms to dry microbursts that 
occur in more tranquil climatic conditions. 
Research in this area continues to gain a better 
understanding of the physics involved- 

The above -eoort also discusses various 
methods r or Jetecting the presence of wind shear 


and alerting the pilot in 3 timely fashion Among 
these 1 s the use of ground based radar3 such as the 
system called NEXPAD. a planned network of Doppler 
radars to be operated by the National Weather 
Service, and a proposed Doppler radar system called 
TDWR that would be placed specifically at principal 
airports. The report also proposes the use of 
airborne Doppler radars, and it is this topic that 
we wish to explore in this paper 

To be more precise, we shall propose a method 
for extracting information from Doppler radar re- 
turns that can be used to determine the presence 
of wind shear. Issues to be discussed will In- 
clude modeling in terms of stochastic differential 
equations, estimation in terms of minimum-mean- 
square-error (MMSE) filters, and detection in 
terms of a first-passage problem. In addition, we 
shall Include some discussion of the hardware 
architecture needed to process the data in a 
timely and efficient manner. 

2. MATHEMATICAL DESCRIPTION 

The method of solution we shall Investigate is 
based on the classical disruption problem In the 
mathematics literature, this is sometimes called 
the first-passage problem or the exit problem. 

More precisely, one monitors the evolution of a 
stochastic process and attempts to determine the 
first time the process reaches a boundary or 
exceeds a threshold. In terms of detecting wind 
shear, the first time the threshold Is exceeded 
will correspond to the onset of a microburst. Our 
formulation will be based on a paper by Davis [4] 
who shows how nonlinear filtering techniques can 
be useful In the detection process. We now 
proceed to describe the radar signal returns, the 
state space model, and the filtering and detection 
algorithm. 

Radar Signal Returns 

The radar signals used In weather applications 
are composites of signals from a very large number 
of scatterers (e.g., hydrometeors) each of which 
can be considered a point target, collectively 
they describe a so-called distributed target. 

After the round trip propagation delay, echoes due 
to one pulse are continuously received over a time 
interval equal to twice the time it takes the 
pulse to propagate across the volume containing 
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the scatterers The echoes can be expressed as 
(see, e g . Van Trees [5] ) 


r(0 - r(2r v /c) 

= j-[FA U e -^ r x/A ]e -jAnr v /A 

22 1 L l l 

t s = 2r v /c is the so-called range 
r„ is the range resolution volume, r. Is 


where 
time,, 

the incremental range within the volume, |Aj |//2" 
is the echo amplitude of the i th scatterer 
located at range r v + r,, and W, is the cor- 
responding range weighting function For a short 
pulse duration, however, this equation does not 
contain any Ooppler shift information because 
weather targets move at relatively 3low speed 
(r, is constant during the time scatterer I is 
illuminated by the pulse). Therefore, In order to 
include Doppler shift Information needed for wind 
shear analysis, one would have to look at returns 
from a number of pulses where each return Is of 
the above form and consists of reflections at a 
number of range gates. By proper sampling, then, 
a sequence of complex video samples r(kT s ) 
separated by the pulse repetition time T s can be 
obtained for each and every range gate This 
sequence would consist of a signal in additive 
noise, I e , 


r(kl 5 j = s k e 


JojdkT., 


V k , k = 0,1, ,M-1 


where w d is the Ooppler shift and M Is the 
number of returns considered at one time. The 
fundamental difference between these equations is 
that the former describes the reflection process 
In space (range time) while the latter describes 
it In time (sample time) at one specific range 
Because our analysis will be based on the sample 
time sequences that contain the Doppler shift 
information, the use of dynamical models that 
Incur abrupt changes are natural candidates for 
wind shear appl ications 


jt;jt e S pace h ade 1 
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The sequence of M points (a) denotes samples 
of the scattered signal at some distance r a in 
the entire volume, (b) denotes samples at dis- 
tance r b , etc Imagine, for example, that each 
return is separated by T 3 = 0 1 ms arising from 
an assumed microburst of 15 km in depth These 
various sequences can be used to generate statis- 
tics for St at the corresponding ranges r v 
within the microburst. For example, one could use 
all M samples in a batch mode to compute the 
radial velocity with a resolution (A/2) (1/MT 5 I, 
where A is the radar carrier wavelength; either 
FFT or covariance methods could be used. 


Now suppose one wishes to estimate the occur- 
rence of the abrupt change in the velocity s^ 
that accompanies a microburst. Consider the 
fol lowing figure 

» r (t> »(*,) 

I ion , 



iom * 



Figure 2 


,n iew of the above discussion let us assume 
that the radar return can be modeled as 


r t - h i s 1 1 ♦ v t 


( 1 ) 


where hls t ) 
A 


3 t 


v t 


“c 


V 


A sin(u c + 4ns t /A)t. 
amplitude i envelope), 
radial velocity (along the radar's 
beam axis). 

Gaussian white noise with zero mean 
and correlation matrix V c , 
carrier frequency, and 
wave I ength 


ke are thus assuming a phase modulation format 
where each measurement r^ represents a return at 
tome fixed range gate, other formats might also be 
appl ’cable More precisely, consider the follow- 
ing figure denoting radar returns from sequential 
pulses of pulse repetition time T s and round 
trip propagation time i s as defined above 


Here we have shown that. In a microburst, the down- 
burst Impacts the ground and the magnitude of the 
velocity v changes abruptly from zero at the 
center of the microburst to some relatively large 
nominal value t v nom , in the absence of a 
microburst, the magnitude and direction of v are 
random 

To formulate a state space model for the above 
situation, let us imagine that before the micro- 
burst. the velocity is a stochastic process x^ 
whose mean Is a relatively small nonzero constant 
and whose variance is a relatively large constant; 
after the microburst, we shall imagine that the 
velocity is a stochastic process b^ whose mean 
is a relatively large constant (corresponding to 
v nom ) and whose variance Is relatively small. 

Thus we can write 


s t = b t z t + x t ( l-z t ) (2) 


where 
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' 0 If 

t < T 

1 lf 

t i T 

5 random 

time c 


(3) 


and T denote 

microburst If we now assume that T Is an exponen- 
tially distributed random variable with parameter 
A. and redefine the observations r t in (I) a 3 
dv t - r t dt. then a possible model might be 


state 


dx 


O' 


dz 


\ 

dt + 

dd 


.0 

dW t 



+ 

dM t 




dB t 


0 

0 



dt 


(4-a) 

(4-b) 

(4-c) 


observation ov^ = A sin(w c * 4nfb^Z(. 

* x t (!-Z t )l/Alt + dV t (5) 


where w t = Brownian motion (mean zero and variance 
W c t), 

M t = martingale associated with a Poisson 
process z^ stopped at its first 
jump time. 

B t = Brownian motion (mean zero and variance 

8 c t). 

V t = Brownian motion (mean zero and variance 
v c t ) . and 

all Brownian motions and initial conditions are 
mutually independent 


The parameters F and G must be chosen to 
achieve the desired statistical properties men- 
tioned above and. at the same time, must guarantee 
stable solutions when their discrete approxima- 
tions are Implemented on a digital computer The 
stability issue will be discussed later in the 
section on simulation Considering the statisti- 
cal properties required of x^, its solution, 
with x Q assumed deterministic, implies the 
following mean and variance - 

tX(- = e'^XQ, var x^ = (W C /2F) ( I — e" j 

Since in reality we are unable to have t » - to 
obtain a large constant variance, we shall choose 
F to be relatively large (and W c even larger). 
However, such a choice drives the mean to zero, 
which is not desired. Therefore, to compensate 
£ or tnls. the x t used in (5) will actually be 
the X(- from (4-a) plus a constant nominal value 
of x nom . On the other hand, since the b^ pro- 
cess has a mean and variance with the same form as 
that for x t , and the desired statistical proper- 
ties of b^ are opposite those of x^, we choose 
G to be relatively small 

Filtering and Detection Algorithm 

Using the state and observation equations in 
(4) and (5)_ above, we can generate MUSE estimates 
x. z, and b of the respective states x, z, and b 
and use them to detect the onset of wind shear. 
More specifically, recalling the definition of 
Z(- in (3). it turns out that 


z t = Prft 2 T|y s . OisSt] (6) 

l.e., the conditional mean z^ Is precisely the 
probability that wind shear has occurred. How- 
ever, because our model represents a system of non- 
linear stochastic differential equations. It is 
not possible to derive the optimal estimate z in 
terms of a f ini te-dlmens lonai realization. There- 
fore, this so-called moment closure problem 
compels us to seek a subopt I ma I filter Once 
specified, then we can use the suboptlmal estimate 
Zt of zt to generate an estimate T of T, the 
first time of occurrence of a microburst. by 
choosing some threshold k e (0,1) and setting 

T = min{t lz t i k) (7) 

It is the form of (7) that causes us to de- 
scribe the problem as one of first-passage Its 
successful use as a detection algorithm will 
depend heavily on the quality of the estimates of 
z provided by the filter, the k parameter is also 
important in that it is related to the performance 
issue and the concomitant false-alarm and miss 
error-probabilities. With this In mind, we have 
opted to compare *"he use of the extended Kalman 
filter (EKF) with the truncated second order 
filter (TSOF) as presented, for example, by 
Jazw inski [6, chap 9] A slight modification in 
the error-covariance equations, however, is 
required because of the martingale M t in (4-b) 
associated with the purely discontinuous process 
Z(- In any case, with overbars denoting 
suboptima I estimates and P denoting the 
error-covariance matrix, we have the following - 



where IN = dy - (h + i • SOT)dt, innovations 
SOT = trace [PH], second-order- term, and 
H = 3^h/3(")^, Hess - an matrix. 


In deriving these equations, use Is made of the 
fact that the definition of z^ In (3) Implies 
z t = z t . Also the term A(l-z) in the 
error-covariance equation is a result of the 
quadratic variation of the martingale M t in 
(4-b). Finally note that if the term SOT is zero, 
then the preceding equations define the EKF - , 
otherwise, they define the TSOF. 

3 SIMULATION 

The approach used in discretizing the con- 
tinuous-time equations given in the previous 
section Is based on Euler's approximation, see, 
for example, Franklin and Powell [7, ch 3] Since 
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the equations In (4) are decoupled, we can Illus- 
trate the approach by discretizing (4-a) to obtain 


4. IMPLEMENTATION 


y k+l = 11 “ Fi)x k + w k A 


where w k Is a zero-mean Gaussian white noise 
sequence with correlation matrix W c To be 
stab'e. one requires that 1 1 — Fi I < 1 . i e., 0 < F 
< 2/4 for a sampling interval of length 4 
which, for us, has been chosen to be 1 Analogous 
statements hold for { 4— c ) Therefore, recalling 
the discussion on parameter selection in the 
previous section on state-space models, we choose 
f = 1 90 and G - 01. 

r ne results obtained thus far from computer 
smuljrions have been encouraging and lead us to 
ccnu'ude Hut the overall approach is a sound one 
«s i s so often the case in such simulations, care 
must be taken in selecting values of state and 
observation noise statistics and initial erroi — 
co/ur lances to prevent filter divergence, and we 
are continuing to investigate adaptive techniques 
to remedy this problem However, in those cases 
where tne filtered estimates did converge, the 
decision that wind shear had occurred was made 
within a few sampling intervals of the actual 
Occurrence It is anticipated that a complete 
pres-ntat i on of these simulations will be 
roi the onn ng in the near future 

we also plan to simulate microbursts through 
the .b» of a model based on the fluid continuity 
equat r on 


7 • (ov) - 0 


with vertical hydrostatic equilibrium Assump- 
tions include i nv i sc id flow with no heat Input 
(dry microburst) An azimuthal ly symmetric solu- 
tion in cylindrical coordinates, with v z and 
v r denoting vertical and radial velocity compo- 
nents, respectively. Is 


v 

z 


be 

1+bH 


-r 2 /2A 2 (e 


-z/b 


-Hz 


e 


) 


2 2 

2A 2 c l - e~ r /2A -z/b 
v r 2 r 

Here, H is the actual scale height, A is the 
initial radius at the top of the column, c is the 
initial velocity or strength, and b is a parameter 
affecting the shape (outflow height). Through a 
s i i gfit modification, rhe model can accommodate 
ssvirnietr i ca I mi crobursts with an essentially 
arbitrary vertical shear profile, the addition of 
a vortex ring to more closely simulate an actual 
microburst, and wet microbursts. Given the 
present state of knowledge concerning wind shear, 
this model produces geometries that are In excel- 
lent agreement with those that have been actually 
observed Futhermore. with this more realistic 
model, we shall be able to simulate radar returns 
-ruin several different range gates 


In implementing the actual system, several 
essential parameters regarding the hardware and 
data acquisition scheme will have to be deter- 
mined. These include Items such as range and 
antenna weighting functions, antenna gain, size of 
the resolution volume scanned by the beam, re- 
ceiver bandwidth, and pulse repetition time. In 
this paper we shall only discuss the latter where, 
ideally, it should be chosen long enough so that 
returns from consecutive pulses do not overlap and 
yet short enough so that samples obtained at one 
range are well correlated for Information extrac- 
tion purposes. We now elaborate. 


Since the system studied here will be used in 
airborne applications, the processing rate, and 
thus the alert time, Isa crucial factor The ‘ 
processing rate has to be fast enough so that the 
pilot can be warned well In advance to avoid wind 
shear. This rate can be Improved by overlapping 
computation time with scanning/sampling time. In 
this context, the system Is thought of as two 
separate units. The first Is a radar that Is 
continuously transmitting and receiving signals 
and a sampler that is sampling the received 
signals and storing them In memory banks according 
to their range as shown in Fig. 3. The second 
unit Is a processor that is computing estimates 
and making decisions While the processor is 
processing returns from the first M pulses, the 
returns of the second M pulses are being sampled 
and stored in memory In order to do this over- 
lapping, the processor must process the data at a 
rate higher or at least equal to the rate at which 


echo amo lit-.de 
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the data are received For example. IF M - 256 
and T s = Ims. then this rate is 1/25 6 kHz. 
guch a high rate may not (depending on the value 
oF N) be achieved with a serial processor 
Instead, one can use a number of Identical 
processors that run in a parallel fashion, see 
Figure 3 Each processor is small, fast, and 
responsible for all the computation required at 
one range gate The computation consists of 
spectral estimation, filtering, and abrupt jump 
monitoring In addition to achieving a high 
processing rate, this scheme is advantageous in 
that N decisions, regarding the presence or 
absence of wind shear, at N range gates will be 
made all at once This Is a very Informative 
aspect because in a wav It reflects the dimension 
of a microburst along the beam axis and It allows 
for rejection of false alarms on part of the indi- 
vidual processors, I e., a microburst should 
trigger the output of more than one processor. 
Thus, a voting algorithm, based on prespecified 
criteria, can be Implemented to produce one final 
dec i s i on 

A large number of range time samples, however, 
imposes a limitation on our processing scheme. A 
large N requires a large number of processors 
wnich could be Impractical due to space and power 
limitations on board an aircraft. In such a case 
the N discrete sequences can be divided into 
groups each of which Is J sequences long, where J 
is a submultiple of N Every group will be pro- 
cessed in a paral lei fashion on J processors, 
while the different groups are fed serially. Even 
in this mi <ed processing, an important improvement 
in the processing rate would be achieved 

5 CO NCLUSI ON 

The preliminary findings described above 
srrcnqlv suggest that the first-passage analysis 
we r ave adopted provides a viable approach to the 
-icb'em of identiryinq the onset or wind shear 
e reva'i's to date center on an idealized 
s "nitric jeonetrv for the wind pattern and have 
ne-i 1 e; r ed noi.e processes specifically related to 
g u"d c'utte’ It is not expected, however, that 
' r- relaxut on or these constraints will alter the 
-a=i_ conclusion future work will incorporate 
eel data and thus will be able to determine the 
u If mate efficacy of the approach 
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APPENDIX B 


Simulated results of state estimates 

Representative state estimates utilizing the model in (4‘) and (5‘) and 
the filter in (10 1 ) and (IT) (with SOT = 0) are shown below for one range 
gate. The filter is implemented in terms of a square-root algorithm; the 
computer program that produced these results on a CDC Cyber 180 computer is 
also given. 

The values of the parameters used in the simulation are as follows: 


A = 1.0 

v c = 1.0 

W c = 0.5 
B c = 0.5 

u c - io ' 5 

X = 0.1 
v = 0.15m 

In each of the four graphs 


G = 10‘ 4 
H = 10" 4 
K : = 32 (s M) 


p u = l.o, p 22 = 1.0 

P 33 = 20.0 

P 44 = °- 5 

P ij = °- 1 » 1 # J 


\ n = 10.0, x n = 10.0 

0 * 0 

b„ = 40.0, b = 30.0 

*0 ‘ *0 = 0 

that follow, the true value of the state is 


denoted with a solid line and its estimate is denoted with a dotted line. 
For example, for state x, the wind velocity before wind shear, the true 
value is Sx and its estimate is Fx; for state z, the true and estimated 
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values are Sz and 
samples, with the 
at time increment 


Fz, respectively; etc. Also, each graph shows 
jump time (indicating the onset of wind shear) 
32. 


63 time 
occurring 
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' . K 1 . I , J , 1 1 
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S X f r : 1 27 ) » S 




n A T A 


k^}(j:l7'),f)jio:i27)OF(0:1Z7),F5(c: 127),P14(0:127), 
P2A(":i27),P34{c:i27)t°4‘iOsi2 7)*SS(0:i2 7), 


(4,4),*s(A f 5),2(4),T(f),V(<0,N(A), 
r ,V'~«\»CK i K ( 1 1 , E R , * X »'<, = E,CF,SA»Fi,F,G,5l' » E C » 

C , f 4 , C t »Pl»32» c ?*R i t»k'5»'t»SI»T, v 'PiAL»l. l $X»L;S&»UFX«UF8 
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R F * r ( 1 * . 1 r ) y l 
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yoirc(?n,i5 ) 


> FPRf*AT(5X) 
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RFAC(10,' , 0)VC,«lC,eC,UC,A,F,G*h»ER,Pll<0),P22(C),P3 3(C>,P44(0) 

cnoMATfn.il 


40 
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epcf'ATf • 


1 2' 



' ( F» ) ,5t2X»'A= , ,c c .:, 
= ’♦ 12///) 


* .5 

CALL PAfOET(S) 
WRITE (20, IPO) 
cpRf'ATf ' f'.'X 

1 • Pll' , c < 


« , F2'.6X. , 53 , ,t<, , F?',fcX, , FSSfcX, 


<♦ ' P22 ' , 5> , * P23* ♦ ! X, ' P4* • ) 


ufx»ip,o*f 

USX*'JFX 
L'FB*3P .0*0 




CF»CF*CO 
00 1000 I*0» 3 


PP 5 00 K = C , 3 1 
TF ( * ,)E. C ) 

r 'ui=pmc 



OP T C 112 



L , (l,3)=(P13(0)-l'(l,4)*t'(3.M*C(M>/C<3) 
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0(1) -P 11 ( C L-1C U » l,(1.2)»»2.C-M? )*J( 1*3) **2. 3- 3( 4) *i; < 1 t f_) **2L. t 
T ( c ) = '■ 











T( 7) 

T( 8) *l'C 
•J ( l » 1 ) - 1 . C 
IJ ( 2 « 2 ) * 1 . C 
U ( 3 * 3 ) * 1 . C 
C ( 4 , 4 ) * 1 . 0 
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10? J*K+I*32 

''2 = 1.0*F2(^) i "x 
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q 2 = c if f ( ) 
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3TBVf Y r "PTT 
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)K=50tT(-?.C*LC«ALCG(t7))SCCS<C0*R8) 
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OF<K+l)-A*SlN{CF*REALCK + l)+SA*SS<K«-in*VK 
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Sl=l*?( K+l) -HK 

S2=((P11(K)*F**2.C+WC)’*C2+EF*F*F12(K)*D3+C*F*P13(K)*04)*C4 


1 

1 

♦F*F+P 1 4 ( K ) *C6 

S3* ( ER *F *P 12 (K >*02*< ER1*( l.O-FZ(K) >♦ P22 ( K ) *E R**2 .0 ) ♦D3+E R* 
G*P23(K )*04)*C4-»ER*H*P24(K)*C6 


1 

$4s(0*c*P13<'<)*C2+S = *G*P23(K)*0 3*(P33(K)*G**2.O3C)*;:4)*C4 
♦ t-*G*P 24 (K ) 7Cb 
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(P44(K)*F**2.0+UC)*Cfc 
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F5(K + 1 )«-FSm*H«-S5/WI*51 



DO 11C J 1* l , 3 

ICS 
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110 
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T ( 1 1 ) =C C 1 1 > 
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CONTINUE 
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CONTINUE 
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C 

CONTINUE 
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CONTINUE 
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CONTINUE 
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